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We present the first fully self-consistent three-dimensional model of a neutron star’s magnetic held, 
generated by electric currents in the star’s crust via the Hall effect. We hnd that the global-scale 
held converges to a dipolar Hall-attractor state, as seen in recent axisymmetric models, but that 
small-scale features in the magnetic held survive even on much longer timescales. These small-scale 
features propagate toward the dipole equator, where the crustal electric currents organize themselves 
into a strong equatorial jet. By calculating the distribution of magnetic stresses in the crust, we 
predict that neutron stars with helds stronger than G can still be subject to starquakes more 
than 10® yr after their formation. 


I. INTRODUCTION 

Neutron stars are of interest not only for the exotic 
states of matter they contain, but also for their magnetic 
fields, which are the strongest in the universe. In the 
case of pulsars, the rotation of the magnetic field pro¬ 
duces beams of non-thermal radiation that can be de¬ 
tected thousands of light-years away. The field can also 
affect the dynamics of the star itself, through the mag¬ 
netic stress it exerts on the star’s solid outer layer, or 
“crust”. In neutron stars with especially strong magnetic 
fields (known as magnetars) the magnetic stress can be 
strong enough to fracture the crust, producing a star- 
quake [I|. This is believed to be the mechanism behind 
the gamma-ray flares and X-ray outbursts detected in 
these objects [3- 

The magnetic fields of pulsars, inferred from spindown 
measurements, can be anywhere in the range 10®-10^® G. 
However, the spindown rate depends only on the large- 
scale component of the field at the magnetic poles, so 
these measurements may underestimate the actual field 
strength in the star. This could explain why some neu¬ 
tron stars produce magnetar-like emissions yet have slow 
spin-down rates There is also evidence that at least 

some neutron stars have stronger magnetic features on 
smaller scales To interpret these observations it is 

necessary to develop a self-consistent, three-dimensional 
(3D) model of neutron star magnetic fields. 

The external field of a neutron star is generated by 
electric currents flowing within its crust and core. The 
ions in the crust form a rigid lattice, and the currents 
there arise purely through the flow of electrons, whose 
dynamics depend primarily on the Hall effect [l0|; this 
situation is commonly referred to as electron magneto¬ 
hydrodynamics (EMHD). Recently there have been nu¬ 
merous studies of EMHD in neutron star crusts [HHil- 
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However, because of difficulties solving the EMHD equa¬ 
tions computationally in 3D, in all of these studies the 
magnetic field was assumed to be axisymmetric. Al¬ 
though there have been 3D studies of EMHD in simpli¬ 
fied periodic-box geometry [20l - l^ , the results cannot be 
used to predict the global-scale field morphology in a real 
neutron star. Currently, the only predictions regarding 
the global-scale field come from axisymmetric models. 

Using one such axisymmetric model, Gourgouliatos 
and Gumming [l8| found that an initially dipolar mag¬ 
netic field evolves towards a quasi-steady configuration 
that they called a “Hall attractor”. After ~ 10® yr the 
electric currents in the crust are concentrated in a narrow 
jet around the dipole equator, producing a strong belt of 
poloidal field. However, it is unknown whether this ax¬ 
isymmetric attractor would be stable to 3D perturbations 
in a real neutron star. 

Here, we present for the first time global 3D numerical 
simulations of the magnetic field in a neutron star crust. 
We use a pseudo-spectral code that allows us to perform 
simulations that are not only fully 3D, but also higher 
resolution than any of those previously presented even in 
2D. We describe how the Hall attractor is modified in 3D, 
and discuss the magnitude and distribution of magnetic 
stresses within the crust. 


II. THE MODEL 

We work in a reference frame that corotates with the 
neutron star’s crust. Because the ions are fixed within 
the crust, in this frame the electric current, J, depends 
only on the electron fluid velocity, v. The current is also 
directly related to the magnetic field, B, via Ampere’s 
Law, and we have (in Gaussian cgs units) 

J = -env = V X B, (1) 

where n is the electron number density, e is the elemen¬ 
tary charge, and c is the speed of light. If the crust is 
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either sufficiently cool, or close to isothermal, then the 
magnetic field is frozen to the electron fluid , and 
evolves according to the equation 


where 77 is the magnetic diffusivity. The two terms on the 
right-hand side of Equation ([2|) represent the Hall effect 
and Ohmic dissipation, respectively. 

In this study we approximate the structure of the crust 
as fixed and spherically symmetric, and so the electron 
density n and magnetic diffusivity rj are fixed functions 
of radius r. We do not therefore take account of any 
deformations in the crust arising from magnetic stresses. 
However, we can use our results to calculate the magnetic 
stress and determine whether it would be large enough 
to induce fractures in the crust. We take the top and 
bottom of the crust to be the spherical surfaces r = R 
and r = 0.9i? respectively, where R = 10km is the radius 
of a typical neutron star. For n(r) and rjir) we adopt 
the same analytical profiles used by Gourgouliatos and 
Gumming [ 2 ^. Specifically, 
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where no = 2.5 x 10^^cm“^ and rjo = 4.0 x 10“'^cm^s“^ 
are the values a.t r = R. These profiles are only rough 
approximations to the (highly uncertain) profiles in real 
neutron stars, but fortunately our results are not sensi¬ 
tive to the specific profiles used. 

Finally, we must impose boundary conditions on the 
magnetic field at the top and bottom boundaries of the 
crust. At the top we impose vacuum boundary condi¬ 
tions, i.e., we match to a current-free field outside the 
star. At the bottom we impose either the idealized 
boundary conditions used by Gourgouliatos and Gum¬ 
ming [^, which are Br = 0 and Jj. = 0, or the more 
realistic boundary conditions used by Hollerbach and 
Rudiger [l^, which model the star’s core as a type-I su¬ 
perconductor. From here on we refer to these two sets 
of conditions as GG and HR respectively. In order to 
implement the HR boundary conditions, it is convenient 
to make a small modification to the electron densi ty p ro- 
file 77 (r), as described by Hollerbach and Rudiger to 
make 1/n vanish at the bottom of the crust. We there¬ 
fore use a slightly modified density profile n{r), defined 
as l/7T,(r) = l/7i(r) — l/7i(0.9i?). 

The relative importance of the Hall effect and Ohmic 
dissipation terms in Equation @ depends on the 
strength of the magnetic field, Bq say, and is quantified 
by the Hall parameter, H = cBo/{'^i^Qno'qo)- In a typical 
magnetar, with Bo = 10^^ G, we have H ~ 50, implying 
that the Hall effect dominates the dynamics of the mag¬ 
netic field. In that case we expect the field to evolve on 
the Hall timescale, 

iHaii = 1.6 Myr, (3) 


where we have assumed that the characteristic length- 
scale for the magnetic field is i? = 10 km. If the field has 
strong, small-scale features, then these will evolve on a 
shorter timescale. 

To solve the EMHD equation ([2]) we have adapted the 
3D MHD code PARODY, developed by Dormy et al. [ 2 ^ 
and Aubert et al. [ 2 ^. The code is pseudo-spectral, and 
uses spherical harmonic expansions in latitude and lon¬ 
gitude, and a discrete grid in radius, making it perfectly 
suited to solving problems in spherical-shell geometry. 
The results that we present here have a resolution of 128 
grid-points in radius, and spherical harmonics up to de¬ 
gree I = 100. We have benchmarked the code against 
previously published axisymmetric results [lj,[^[^ and 
find excellent agreement in all cases. 


III. RESULTS 


A. Robustness of the Hall attractor 

To determine whether the “Hall attractor” seen in ear¬ 
lier axisymmetric simulations is robust against 3D per¬ 
turbations, we have repeated one simulation of Gourgou¬ 
liatos and Gumming (23|, using the same boundary con¬ 
ditions and dipolar “Ohmic eigenmode” initial condition 
for the axisymmetric component of the magnetic field. To 
this initial axisymmetric poloidal field we add a low am¬ 
plitude, small-scale 3D perturbation with both poloidal 
and toroidal components, comprising spherical harmon¬ 
ics of degree 20 ^ ^ 40. As shown in Figure [U the mag¬ 
netic field evolves towards a state broadly similar to the 
axisymmetric Hall attractor. However, significant non- 
axisymmetric features persist in the simulation even on 
long timescales (comparable to the global Hall timescale 
^Haii)- These features would rapidly decay by the action 
of Ohmic dissipation alone, so their longevity can be at¬ 
tributed to the Hall effect. In fact it has previously been 
shown that, in the presence of a strong density gradient, 
the Hall effect can sustain or even amplify small-scale fea¬ 
tures in the magnetic field [l^. The length-scale of the 
longest-lived non-axisymmetric features is comparable to 
the largest scales present in the initial 3D perturbations 
(degree I = 20), suggesting that they form by an up-scale 
transfer of magnetic energy. 

The convergence of the global-scale magnetic field to 
the attractor state occurs on the Hall timescale ([3]), after 
which the evolution proceeds on the slower timescale of 
Ohmic diffusion. This is in a gree ment with results from 
axisymmetric simulations [ij, [l^ . (In the simulations of 
Marchant et al. [I^, by contrast, the field was initialized 
already in a Hall equilibrium state, and therefore evolved 
on the Ohmic timescale throughout.) The poloidal mag¬ 
netic flux becomes concentrated into a “belt” around the 
equator of the magnetic dipole, implying a strong equa¬ 
torial jet of electrons within the crust. The angular ve¬ 
locity in this jet is approximately constant along each 
poloidal field line (see Figured]). We do not find that the 
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FIG. 1. (a) Lines of the azimuthally averaged poloidal field and electron angular velocity at successive times, illustrating 
convergence to the Hall attractor. The dashed lines indicate the boundaries of the crust, (b) The radial component of the 
magnetic field at the surface at early and late times in the same simulation. 


Hall effect enhances the diffusion of the magnetic field, y 
originally suggested by Goldreich and Reisenegger [lOl| . 
However, by concentrating the poloidal flux around the 
equator, the Hall effect does lead to a more rapid decline 
in the strength of the radial field at the poles. 

The magnetic field is primarily poloidal, with a weaker 
toroidal component in the inner crust. However, the GC 
boundary conditions used in this simulation impose zero 
toroidal field both at the top and bottom of the crust. 
In a real neutron star, currents can flow between the 
crust and core, generating significant toroidal fields. We 
have therefore repeated this simulation using the more 
realistic HR boundary conditions, which allow for finite 
toroidal field at the lower boundary. Figure [2] compares 
the azimuthally averaged magnetic field at the same time 
in both simulations, revealing that although the mor¬ 
phology of the poloidal field is modified, a similar Hall 
attractor state still exists. The toroidal field in the sec¬ 
ond simulation is much stronger, as expected, but re¬ 
mains weaker than the poloidal field. This toroidal field 
is associated with a meridional flow of electrons that is 
equatorward at the star’s surface, and drags the poloidal 
magnetic field lines. As a result, the equatorial belt of 
poloidal flux, and the corresponding electron jet, is even 
stronger in the simulation with HR boundary conditions, 
and is pushed slightly deeper into the crust. Figure [3] 
illustrates the development of the equatorial belt, which 
manifests as bands of positive and negative Br on either 
side of the equator at the star’s surface. 


B. The effect of field strength 




FIG. 2. Poloidal fieldlines and contours of toroidal field at 
t = 1.2 Myr. The left and right panels show the cases with 
GC and HR boundary conditions, respectively. 


surface. The Hall effect may therefore be even more sig¬ 
nificant in neutron stars than previously thought, and so 
we have repeated the same simulation shown in Figure [3] 
with a stronger initial magnetic field (i.e., a larger Hall 
parameter H). Because the evolution of the global-scale 
magnetic field occurs on the Hall timescale (|3]), we expect 
that a stronger field will converge more rapidly to the 
Hall attractor, before the small-scale magnetic features 
have been dissipated by resistivity. This is confirmed in 
Figure SI which shows the surface radial field in three 
simulations with increasing magnetic field strengths. In 
each case the field is plotted at time t = O.fitHaii, which 
corresponds to t ~ I, 0.5, and 0.25Myr respectively, and 
an equatorial jet of electrons has already formed. In the 
case with the strongest magnetic field, the jet is broader 
and more spatially disordered; at later times in the same 
simulation, the jet becomes increasingly laminar, but re¬ 
mains broader than in the other simulations. 


Spin-down measurements of neutron stars tell us only 
the strength of the large-scale magnetic field near the 
poles. However, it is clear from Figure [3] that the field 
strength elsewhere in the crust can greatly exceed this 
observed value. At the time plotted in Figure for ex¬ 
ample, the strength of the surface magnetic field varies 
between ~ 5 x 10^^ G at the dipole axis and ~ 3 x 10^^ G 
near the equator. Near the bottom of the crust, the field 
is stronger by a further order of magnitude than at the 


C. The crustal magnetic stress 

If the magnetic shear stress within the crust exceeds 
the breaking stress of the ionic lattice then it will induce 
a crustal fracture. Molecular dynamics models indicate 
that the breaking stress is approximately 5 x 10“^ of 
the electron pressure, Pe [1^. At each point within the 
crust, the strongest magnetic shear stress is exerted on 
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FIG. 3. Evolution of the magnetic field in the simulation with HR boundary conditions. Top row: at r = R. Bottom row: 

Bcf, at r = 0.9i?. Colorbars are adjusted to the maximum values in each plot, and use a linear scale. 
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FIG. 4. The surface radial field in three simulations with Bo = 10^"^ G, 2 x 10^"* G, 4 x 10^"* G (i.e., H = 50,100, 200). 


a surface that makes an angle of 45° to the local mag¬ 
netic field direction, and this stress is exactly equal to 
the magnetic pressure. Pm- The condition for fracturing 
is therefore Pm ^ 5 x 10“^Pe, where Pm = |Bp/(87r) 
and Pe ~ (37r^n)^/^/ic/(127r^). Fractures are most likely 
near the surface of the crust, where the density and pres¬ 
sure are lowest, and the energy released in a near-surface 
fracture can directly power flares and outbursts. In Fig¬ 
ure [5] we plot the ratio of magnetic pressure to breaking 
stress at the surface of the crust in the same simulation 
shown in the last panel of Figure|4j revealing that patches 
around the dipole equator are susceptible to crustal frac¬ 
turing. These patches are localized in both latitude and 
longitude, as a consequence of the three-dimensionality 
of the magnetic field. They persist on long timescales, of 
order tnaii, and propagate azimuthally in the direction of 
the electron jet. 


IV. DISCUSSION 

Our results indicate that neutron starquakes are most 
common in the vicinity of the electron jet within the 
crust, which typically forms a ring around the dipole axis. 
This jet forms on the Hall timescale, which is dependent 
on the overall strength of the magnetic field. At earlier 
times the structure of the magnetic field is dependent on 
the initial conditions for the proto-neutron star, which 
unfortunately are not well known. The surface field is 
strongest in localized patches around the magnetic equa¬ 
tor, where the local field strength typically exceeds that 
at the poles by an order of magnitude. 

In our simulations, surface magnetic features migrate 



FIG. 5. Lines of the magnetic field generated by currents 
in the crust. The coloring indicates the ratio of magnetic 
pressure to breaking stress at the surface of the crust. A 
portion of the surface has been cut away to show fieldlines 
inside the crust, whose lower boundary is the gray sphere. 


equatorward as a result of the toroidal magnetic field in 
the crust and the associated meridional flow of electrons. 
This migration could potentially be reversed if a strong 
toroidal field of the opposite sign were present, corre- 
sponding to a poleward flow of electrons at the star’s 
surface [^. There is no obvious process that can gen¬ 
erate such a strong toroidal field within the crust, but 
it could be a remnant from the star’s formation l3C| or 
the result of toroidal flux expulsion from the core [3l|. 
To fully describe the processes that can impart a strong 




























5 


toroidal field to the crust, it will be necessary to couple 
our crust model to a model of the superconducting core. 

The flow of heat from a cooling neutron star is inhib¬ 
ited across magnetic field lines, so the strong concentra¬ 
tion of poloidal magnetic flux in the equatorial belt could 
trap heat within the electron jet [33 ■ The breaking stress 
is very sensitive to temperature [28j| . and so the equato¬ 
rial region of the crust may be even more susceptible to 
fracturing than we have found here. We are currently 
extending our model to describe the thermal-magnetic 
evolution of the crust, and its coupling to the core. 
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